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We propose a new mechanism to increase the reactions rates in multistable autocatalytic systems. 
The mechanism is based upon the possibility for the enhancement of the response of the system 
due to the cooperative behavior between the noise and an external periodic modulation. In order 
to illustrate this feature we compute the reaction velocities for the particular case of the Sel'Kov 
model, showing that they increase significantly when the periodic modulation is introduced. This 

OS ■ behavior originates from the existence of a minimum in the mean first passage time, one of the 

OS ■ signatures of stochastic resonance. 
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I. INTRODUCTION 



In the field of chemical kinetics and in particular in the study of autocatalytic systems the enhancement of the 
reaction velocities is a mattep.of primary interest. Some mechanisms, among which we can quote the chemical reactions 
CN ' ft s H s ^ ec ^ ky a magnetic fieldEl, have been proposed to this end. In this paper we propose stochastic resonance (SR) 
( \ an as a new mechanism to increase the reaction rates and diffusion in autocatalytic systems. This phenomenon 
t— I ■ consists in the enhancement of the response of a system by the noise when it is under the influence of a periodic 
input. Thus, it constitutes a cooperative phenomenon between noise and an. oscillating signal, being proper to 
nonlinear dynamical systems, and also to certain class of nondynamical systems! Several magnitudes enables one to 
characterize the presence of SR. The most common of them is the signal-to- noise ratio (SNR), usually used in signal 
detection problems, which measures the response of the system to the input. The existence of SR is revealed by th 
presence of a maximum in the SNR for a non-zero noise lcvcld. However, for our purposes we will pay attention on 
another signature of the presence of SR in the system: the existence of a minimum in the mean first passage time 
"-^J ■ (MFPT) as a function of the input frequencyO. This quantity is, essentially, the relaxation time of the system. Thus, 
£h ! a diminution in the MFPT implies an increase of the reaction rates. _ _ 

Although some works have previously been devoted to show the existence of SR in chemical systemsE3~lij, the 
J-^ ' possibility of improving the reaction rates is an aspect which has not been treated up to know. When we are 
£> . concerned with catalytic systems, the fact that the MFPT, which can be interperted as the relaxation time, reaches 
J^j ' a minimum for some non-zero noise level leads to an enhancement in the time of reaction. ■— ■ 
Along this paper we will focus our attention on a particular example of autocatalytic model: the Sel'Kov modellia. 
To be precise, we will use an implementation of-this model proposed by Ross et al.. To achieve this purpose, we will 
apply to this model a theory we have developed!^ with the objective of deriving analytical expressions for the reaction 
rates and more generally to analyze the stochastic properties of a certain class of nonpotential systems, to which the 
Sel'Kov model beloELgs-j-Making use of this theory we can extend some previous analysis performed in the framework 
of potential systemsO tL3 to the nonpotential case thus, unifying the treatment for both. 

The approach to the problem that we are going to adopted is based on the assumption of weak noise limit and 
on the particular structure of the phase space of the dynamical system. The phase space must be two dimensional 
with three aligned fixed points: one unstable fixed fixed point between two stable fixed point. Moreover, one of the 
variables has to relax faster than the other, in such a way that we can restrict the dynamics, in the adiabatic limit, to 
the line which contains the fixed points. The reduction to an one dimensional dynamics, which is obviously potential, 
allows us to treat the system in the same way as if it were potential. 

We have organized the paper in the following way. In Section II we formulate the kinetic equations for the Sel'Kov 
model and study the behavior of the reaction rates under a periodic external forcing. Section III is devoted to 
discuss the stochastic properties of this model. We show the existence of SR by computing the signal-to-noise ratio 
in Section III, and the distribution of residence times and the MFPT in Section IV. Finally, in Section V we discuss 
the implications of our results in chemical kinetics. 
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II. THE SEL'KOV MODEL 

The equations of the Scl'Kov model, as proposed by Ross et aL0 are given by 

dp x 



dt 



hp A + hPy ~ k 2Px ~ hPxPy (!) 



—37- = hPB + hPxPy ~ hPY - k 4 p Y . (2) 
at 

These equations correspond to a reaction occurring isothermally in a region of constant volume linked to two infinite 
reservoirs, containing species A and B, by membranes selectively permeable. This region contains intermediate species 
X and Y and the catalysts. The two reservoirs are maintained at constant pressures pa and pb whereas in the 
intermediate region the pressure is not constant due to the reaction. Eqs. ([!]) and (Q) give us the evolution of the 
vapour pressures for the intermediate species Px and Py . The external periodic forcing will be introduced through 
periodic variations of pa or ps- 



A. Kinetic Equations 

Let us briefly discuss the method we have used to obtain the kinetic equationsB. The class of two-dimensional 
nonpotential noisy dynamical systems to which our method applies is characterized by the peculiar topology of the 
phase space. For the case concerning us, the dynamics is characterized by the presence of three aligned fixed points: 
an unstable fixed points between two stable fixed points. Moreover, it is necessary that one of the variables relaxes 
faster than the other. In our case this fast variable is_»y whose relaxation time is, approximately, fc^ 1 which is smaller 
than feg , the corresponding relaxation time for pxliE ■— ■ 

Consider a general two-dimensional noisy dynamical systcmEd 

dx ' —* 

— = u(x) + g.£(t), (3) 

where x is the vector whose components are state variables, the field u(x) is the drift, which is suposed to fulfill 
the aforementioned conditions, g is the noise matrix and is a gaussian white noise of zero mean and correlation 
function given by 

= 2D6 t3 6(t-t'), (4) 

with D accounting for the noise level. 

The corresponding Fokker-Planck equation involving the probability density, p(x, t), is 

^ = V-(-u(g)p + V-(Dp)), (5) 

where D = Dg ■ g is the diffusion tensor, with T indicating transposition. 

Let us now assume that the system is potential. In such a case it is possible to write the Fokker-Planck equation 
as a continuity equation 

|-v./, (0, 

where J is the diffusion current given by 

J = -De- U /°V -e// D . (7) 
To obtain this expression we have defined U and /2 ( a generalized chemical potential) as follows 

U = ~f u-dx, (8) 
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9pe U ' U - (9) 



For a potential system the function U denned through (^) coincides with its potential energy. In the nonpotential 
case, however, the value of U will depend on the path of integration we choose and, in general, we cannot achieve eqs. 
@ and ©. 

Let us now assume the adiabatic hypothesis together with the weak noise limit and analyze the consequences that 
these assumptions have on the general properties of the quasi-stationary probability distribution. If we allow the 
system to evolve for sufficiently long time, the probability distribution will acquire two high maxima around the 
stable fixed points (SFP) of the (deterministic) systemn3. In the weak noise limit, the probability distribution will 
be very narrow around the line on which— the fixed points lie. Therefore, we can assume that the fluctuations are 
concentrated along this line. As shown ir£3, the potential for this system is given by the function U taking precisely 
this line as integration path. 

From eqs. (|l|) and (|2|) by equating their right hand side to zero and adding the resulting equations, one can see 
that the fixed points have to lie on the line 

Py = -jrPx H r • (10) 



Having in mind that py is the fast variable, we can introduce ( jlOD in (|l|) to obtain a one dimensional drift, which, when 
equated to zero and for the proper values of the parameters, yields a third order equation with three real different 
roots, corresponding to the steady states of the system. A simple computation allows us to see that eq. (||) when 
integrated along the line ( |l0| ) gives the potential for this drift. 

Before proceeding further, it is convenient to rescale the resulting one dimensional system. To this purpose we will 
introduce the dimensionless variables 

x = — , 11 

PA 



t = k 2 t. (12) 

The corresponding Langevin equation describing the dynamics of the fluctuations along the line defined through 
eq. ( [l0|) is, in nondimensional variables, 

x = a + b - (1 - c)x + dx 2 - ex 3 + f (r), (13) 

where the coefficients are related to the original parameters by 

fci 

a = T~! 

k 2 

b = K 3 ^, 
PA 

d = K(p A —— + 2pA — ) , 

2 k 2 , kik 2 . 

where K = fclPA + fcePB . The potential as a function of this nondimensional variable is obviously given by 

U = -x^ - ^-x 3 + ^x 2 - (a + b)x. (15) 
4 3 2 

In Figs, la and lb we show the fixed points of the Sel'Kov system and the form of the function U, respectively. 
From Kiff. la it is easy to infer that the unstable fixed point (UFP) is Fq, while F± are the SFP (of the 2-dimensional 
model)E3. Fig. 2 shows that the SFP of the total systems are in this reduced version the minima of the potential 
corresponding to the maxima of the probability distribution function, and the UFP is a maximum of U. 

Our next step will be to discretize the dynamics on the null cline . In particular we will obtain theikinetic equations. 
To this end, we define the populations n + ( n_ ) as the population on the right ( left ) of the UFPE3 
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p(x,t)dx (16) 
s[+] 



p(x,t)dx (17) 

SH 

where S[+] (S[— ]) is the portion of the phase space on the right( left ) of a line passing through the UFP and 
orthogonal to the line which contains the fixed points. 

In the adiabatic limit, we can assume that the population is strongly concentrated in a small region around the 
SFP, as suggested by the picture of the probability density we have profiled when the maxima in this long time limit 
is very high. Therefore, in this limit we can imagine that the system reaches a quasi-stationary state in which a 
quasi-stationary diffusion current is established. This current is assumed to be uniform between the two maxima of 
the probability density and, in the weak noise limit concentrates around the line joining the three fixed points 

J(x,t) = J(x,t)S(y - y (x)) = J(t)S(y - y (x))(9(x - x+) - 6{x - (18) 

where yo{x) refers to the line ( |l0| ) written in nondimensional variables and x+ ( x_ ) is the coordinate of the fixed 
point on the right ( left ) of the UFP. 
The kinetic equation for n + is given by 

dn+ f %x = -f V-Jdx. (19) 



dt Js[+] dt j 3 [. t 

By using the divergence theorem and the assumptions about the form of the diffusion current, we have 

= J + Ji(x,t)S(y)dy = JiOM), (20) 

and, proceeding in the same way as for the case of n_, we obtain 

dn- , , . 

= -Ji(x,t) (21) 



dt 



In addition, due to the height of the maxima in the probability density and the weakness of the noise, we can also 
consider that equilibrium in each side of the UFP is reached independently. Thus, the generalized chemical potential 
is given by 

H(x, t) = {fi(x + ,t)9(x -x) + n(x-,t)6(x - x )}5(y - y (x)), (22) 

where Xq is the coordinate of the UFP. Notice that the tensorial character of the generalized potential has been 
removed due to the reduction to a one dimensional dynamics. Moreover, we have defined /i = fin. By using equation 
( p2| ) in ^ we then obtain 

*(£, t) = {^ +e -^ u ~ u +^ D 9(x -x)+ ^-e- { - u - u -y D 6{x - x )}S(y), (23) 

where U corresponds to the integral along the adequate path (|l0|), U+ and C/_ are its values at the SFP ( its minima 
), * = g n p and *± = gnp(x±,t). j-. 

In order to obtain the expression for the quasi-stationary current J\ (t) we follow the same procedure as inll3. Finally, 
we obtain the kinetic equations for n± 

dn + dn_ , , 

-w—nr-*-"-**"*' (24) 

where the kinetic coefficients are given by 

K T = U\u;\u';)V*e-V°- u T)/n. (25) 
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B. Periodically Modulated Sel'kov Model 



Let us now introduce an external periodic forcing in this model by considering a weak periodic perturbation in ps 
of the form 

p B = Po(l + e(t)), (26) 

where e(t) = eosin(oJo^); with eo being a small parameter. The pressure pa is mantained constant. In this case the 
coefficients a and e remain constant while b, c and d become periodic functions, 

b = b Q (l + e(t)), 

c = Co(l + e(i)), (27) 
d = do(l + c(*)), 

Since U, its derivatives and-thc positions of the fixed points are also periodic functions, the reaction rates K± oscillate, 
too. For these we achieveOu, up to second order in the parameter (fro = ^ and under the assumption of weak noise, 

K±(t) = a^ + af 4>o sin ui^t + aftfro sin 2 u>ot, (28) 
where a , af and af are given by 

4 = ~^f^^ D , (29) 

Z7T 

af = - C -^tfe^ D { m - V± ), (30) 

at = l^^tie^ino-n^. (31) 

Here /q is the zero order contribution to (IC/g'j/7^) 1 / 2 . Notice that the existence of a periodic forcing is responsible 
for the time dependence of the reaction rates. 

In Fig. 2a we have plotted K + as a function of time for a fixed noise level. As we can see from this figure, there 
appear periodic oscillations in the reaction rates. This behavior is a consequence of the form of the distribution of 
residence times ( see Section IV). An interesting phenomenon occurs when the noise level is very low. In this case, 
the reaction velocities become very narrow around the maxima and are practically zero at any other time. 

Let us now focus our attention on the time averages of the reaction rates. These are given by 



(K ± ) t = ^-Jj l dtK ± (t) = fte^Xl + ^lim-^Y) 



(32) 



It is important to note that, although in some time intervals the reaction rates are smaller than in the non perturbed 
case, the average over a period of the external forcing is greater than the value of the reaction rate in the non perturbed 
case, corresponding to the same noise level. This can be seen by comparison of Figs. 2b and Fig 2c. 

III. POWER SPECTRUM 

We are now in a position to obtain the stochastic properties of the Sel'Kov model. In particular, we will compute 
the power spectrum and from it the signal-to- noise ratio (SNR) . 

In order to compute these (Quantities, we assume that, in the limit of weak noise, the probability density ( in one 



dimension ) can be written asi 



p{v,t) = n + (t)8 V:V+ + n-(t)8 v<v _, (33) 



where n+ and n_ evolve in time according to the kinetic equations (|24|). The formal solution of these equations is 
found to be 
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n±(t) = g-Ht)Mto)g{to) + I K T (t')g(t')dt'), (34) 



to 



with 



g(t) = J (K + + K_)dt'. (35) 



Introducing the expansion (|2S]) for the reaction rates in ( |34|) and (|3^) one obtains, up to first order in 

ai cos(uj t + $) ai cos(cj t + *) 

w («o+ w o) 7 "o ("o + ^o) 7 

ai sin(uJot + $) 
^(ag+a; 2 )V2' 



(36) 



where $ = arci<7(a;o/wo). The quantity n±(t\vt ,to) is the conditional probability that v(t) is in the + state at time 
t, given that the state at time to w &s vt . From this equation it is possible to compute the-,statistical properties of the 
process v(i). Of particular interest to our purposes is to find its autocorrelation functional 

(v(t)v(t + T )) = lim t ^-oo(v(t)v(t + T)\v to ,t ). (37) 

The conditional correlation function is given by 

(v(t)v(t + r)\v ta ,t Q ) = v + (t + T)v + (t)n + (t + T\v + ,t)n + (t\v to ,t ) + 
v + (t + r)w_(<)n + (< + T\v-,t)n-(t\vto,to) + 
V-(t + r)v + (t)n_(t + T\v + ,t)n + (t\v to ,to) + 

v_ (t + t>_ (t)n_ (t + t \v+, (t\v to , t ). (38) 

Let us make some considerations which will allows us to simplify the computation of the autocorrelation function. 
It is clear from equation (38) that the Fourier transform of the autocorrelation function will depend on t as-well as 
on the frequency. This dependency is avoided by taking its average over the period of the external forcinga. The 
autocorrelation function is then computed up to second order in the parameter <po ~ _D _1 .Thus, in the limit of the 
weak noise D~ x it can be neglected when compared with D~ 2 . Therefore, the only contribution of the first order 
term of v± to the autocorrelation function comes from its product with the zero order term of the product of n's. 
But on performing the average this term vanishes. Moreover, the position of F— is very near to zero. In fact, it is an 
order of magnitude smaller than eo and therefore it can be neglected. Finally one arrives at 



(v(t)v(t + T)\v to ,t ) = r) + n + (t + T\v + ,t)n + (t\vt ,t ), (39) 



where the bar indicates average over t. From equations (38) and ( |39[ ) taking the average and the limit to — ► — oo we 
finally obtain the following expression for the autocorrelation function 

M*M* + r)> = Vl(^) 2 + e— - 2£) + 

Vl^de-^i^^ - M^) 2 " («i) 2 - (-) 2 )(^^) + 
a^a^ai a\ 2 cos(uj t) 

{ «0l ) 2" + 

(— — - + 1 )sm{ujQT)) + 
a Q luq a 

cosM _ a- oa - l a 1 + I + 1 a, 2 

ag +Wp Wp a 2 2 uj q 

Whith these results we can now compute the average of the power spectrum defined as 



G 



W ° / OCO i\Ji _ / 7J7w7 . -,ii T , 



5(fi) = ^ / S(n,t)d* = / (w(i)w(t + T))e- liiT dr, (41) 
2t Jo J-og 

where the last equality follows from the commutative character of the average and Fourier transform. After Fourier 
transforming (40), we obtain 

W) = V 2 AS(SI) + 2^-^(1 - + 2.^4(^1 - 

k^mn - c„). (42) 

"o 2 a Q 2 ll> 

In this expression, the fraction of the total power ua the broadband noisy part of the spectrum, which usually is a 
small fraction of the total power, has been neglectedu. Note that the power spectrum contains a term proportional to 
5(H). The appearance of this term is due to the asymmetry of the potential that originates a mean probability current 
between the two stable states, i.e. towards the deepest well of the potential. This corresponds to the (deterministic) 
approach to the equilibrium state. 

From equation (42) the signal to noise ratio, SNR, can be obtained as a function of the noise level D, by making 

fl = UJQ. 

£aJ(3 a - (S i)2 _ + i (££)!) 
SNR = °- ^ — . (43) 

I a o 

a a 

This quantity has been plotted in Fig. 3 as a function of the noise level or input noise. It exhibits a maximum at 
a certain value of D, thus indicating the existence of stochastic resonance. 



IV. ESCAPE-TIME DISTRIBUTION AND MEAN FIRST PASSAGE TIME 

In order to compute the escape-time distribution ( ETD ) around the fixed points of the system we will assume the 
existence of an absorbing barrier between thema. This implies that the kinetic equations reduce to 

^ = -K ±n± . (44) 
dt 

In this case, and using initial conditions n±(t = 0) = 1, the solution of eq. ( f44| ) can be written as 

n±(t) = exp{ / dzK±(z)}. (45) 

^0 Jo 

On the other hand, the ETD around the fixed point is given by 

p± = (46) 

Therefore, the ETD reads 

P± = A±cxp{ / dzK±{z)}. (47) 

Jo 

Performing an expansion up to second order in the parameter (f>o we obtain 

P± = i a o + a i 4>o sm ^ot + a 2 00 sm2 ^ot) x 

exp{-aftt + -L0o(l - cos^) + — 0o( T 1 ( 48 ) 

From this equation we can see that a coherent response of the system is generated by the input signal in the sense 
that the distribution is modulated by the signal. We have plotted this quantity in Fig. 4. 
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As follows from the definition of the ETD, the quantity p±dt corresponds to the probability that the absorbing 
boundary to be reached when the evolution of the system comes from the fixed point F± . The MFPT of leaving F + 
is given by 

/>oo />oo 

(T) = / dttp+{t) = I dtn + (t). (49) 
Jo Jo 

where use has been made of eq. (fl6|) and an integration by parts has been performed. Using eq. ([is]) with eq. (|2S| ) 
we obtain 

(T) = I" dtexp{- — (a+t + a+Ml-oosLJ t) + a+ ( l>l( ( ^- S -^^)}. (50) 
Jo ^0 Z i 

Taking into account that eo is a small parameter, we can perform a Taylor expansion of the exponential. Up to second 
order we obtain 

-; — r = 1 — OlT 0n i k 

(To) 1% (a+)2+^ 

2 l lM ((^) 2 +4c 2 )((a+) 2 +^) 

1 a+^l 4w 2 

2 a + K) 2 + ^- 15 > 

where (To) is the MFPT in the absence of external forcing. We have plotted this quantity in Fig. 5. 



V. DISCUSSION 



In this paper we have proposed a new mechanism to increase the reaction rates or equivalcntly to reduce the reaction 
times in multistable autocatalytic systems. This mechanism is based upon the existence of SR in this kind of systems, 
which consists essentially in the enhancement of the response of the system to a periodic input and the diminution of 
the MFPT. 

We have found that the introduction of a weak periodic forcing improves the reaction velocities. As concluded from 
Fig. 2b, 2c and 3 the values of the input noise for which this effect is more important correspond, roughly, to the 
interval of noise intensity where SR occurs. To be precise, for D = 0.19, corresponding to the appearance of SR in 
the system, there is a significant increase in the reaction rate of about 40 per cent. 

Another interesting result about the behavior of the reaction rates concerns its temporal dependence, which can 
be explained on the light of the behavior of the ETD. Comparing Figs. 2a and 4, one can observe that the reaction 
rate reaches its maximum value when the corresponding ETD is at a maximum and vice versa, hence, when this 
probability reaches a maximum value the corresponding transition rate must also be maximum. 

With regards to the time dependence of the reaction velocities, we have found that, in the limit of very weak noise, 
they become very narrow around the maxima and near zero for any other time. In order to explain this fact let us come 
back to the picture of the quasi-stationary probability density, solution of the Fokker-Planck equation, described in 
section II. The width of the oscillations around the maxima in the reaction rates, as well as in the ETD, are due to the 
transitions induced by the noise. In a situation of very low noise intensity these transitions are practically absent. In 
this case the transitions between the fixed points have their origin, mainly, in the deterministic dynamics which allows 
the system to go from one fixed point to another, as reflected in (42) by the presence of the term proportional to <5(f2), 
and in the periodic modulation of the probability density. This makes the height of the maxima of the probability 
density to change along the time. Consequently, the higher of them becomes shorter than the other one and this 
produces a current to the new most probably position. This last effect is responsible for the periodicity observed in 
the transition rates, whereas the transitions between peaks originate from the peculiar deterministic dynamics. 

It is important to note that these features about the behavior of the probability density are reproduced by the 
function U. For instance, the aforementioned periodicity of the transition rates can also be understood in terms of 
the periodic modulation of the function U, which plays the role of a potential. 

On the other hand, due to the existence of SR in the system, the MFPT, which can be resembled to the reaction 
time, exhibits a minimum as a function of the frequency of the periodic input. Moreover, in the context of transport 
phenomena a MFPT can be understood as a relaxation time which, in essence, is the inverse of a diffusion coefficient. 
Therefore, a minimum in MFPT implies an improvement in the diffusion of the system. As one can see from Fig. 5 
the diminution in the MFPT with respect to the case in which the system is not under the action of a periodic forcing 
is significantly about 20-25 per cent. 
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FIGURE CAPTIONS 

Figure 1.- 

— (a) Fixed points of the Sel'Kov model for a = b = 0.5, c = 1.1, d = 3 and < 

— (b) Asymmetric one dimensional potential as a function of x. 

Figure 2- 

— (a) K + as a function of time for D = 0.19 and eo = 0.03. 

— (b) K + in absence of periodic modulation as a function of D. 

— (c) Temporal average of K + (t) as a function of D for eo = 0.03. 

Figure 3.- Signal-to-Noise Ratio as a function of the noise level D for e = 0.03. 
Figure 4.- Escape time distribution as function of oj^t for D — 0.19 and e = 0.03 
Figure 5.- Mean first passage time as a function of luq for D = 0.19 and e = 0.03 
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Fig. 3 




Fig. 4 




Fig. 5 




